Geometrical Frustration: A Study of Ad Hard Spheres 
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The smallest maximum kissing-number Voronoi polyhedron of 3d spheres is the icosahedron and 
the tetrahedron is the smallest volume that can show up in Delaunay tessalation. No periodic 
lattice is consistent with either and hence these dense packings are geometrically frustrated. Be- 
cause icosahedra can be assembled from almost perfect tetrahedra, the terms "icosahedral" and 
"polytetrahedral" packing are often used interchangeably, which leaves the true origin of geometric 
frustration unclear. Here we report a computational study of freezing of Ad hard spheres, where 
the densest Voronoi cluster is compatible with the symmetry of the densest crystal, while polyte- 
trahedral order is not. We observe that, under otherwise comparable conditions, crystal nucleation 
in Ad is less facile than in 3d. This suggest that it is the geometrical frustration of polytetrahedral 
structures that inhibits crystallization. 

PACS numbers: 64.70.Dv, 82.30.Nr, 81.10.Aj, 64.70.Pf 
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Most glasses form under conditions where the ther- 
modynamically stable state of the system is crystalline. 
Good glass formers should therefore be poor crystalliz- 
ers. Geometrical frustration is one of the factors that 
may prevent the formation of the ordered phase and 
therefore help physical glass formation There is also 
evidence that such frustration increases the height of 
the crystallization- nucleation barrier of liquid metals • 
Isotropic simple liquids are often considered frustrated 
because the five-fold symmetry of the liquid icosahedron 
cannot pack as a regular lattice. This scenario contrasts 
with what happens in a fluid of 2d disks, where hexago- 
nal order is both locally and globally preferred and where 
crystallization is particularly easy. 

Several physical mechanisms have been proposed to 
support the formation of icosahedra. On the one hand, 
Frank, considering the optimal way for kissing spheres to 
cluster around a central one, found the icosahedron to 
be more stable than the cubic lattice unit cells for the 
Lennard- Jones model Q . Though the original argument 
relies on the energetics of spurious surface effects [4j, 
mean-field studies correcting for solvation leave the result 
unchanged [!, 0|. The icosahedron, the smallest max- 
imum kissing-number Voronoi polyhedron, is optimally 
packed. It offers the most free volume to surface spheres, 
so it is also preferred entropically. On the other hand, the 
polytetrahedral scenario ascribes the presence of icosa- 
hedra to their facile assembly from quasi-perfect tetrahe- 
dra, themselves the smallest Delaunay decomposition of 
space 7, 8] . But is it the packing of Voronoi polyhedra or 
the packingof Delaunay hyper-triangles that counts? Ex- 
periments jlidclllll and simulations [H|,[l3[ only manage 
to identify icosahedral order in limited quantities, even in 
deeply supercooled systems. Recent studies indicate that 
liquid polytetrahedral order is a lot more varied 3, HI 



than the icosahedral picture suggests. Yet, because of 
the geometrical ambiguity, the equation of the icosahe- 
dron with frustration is difficult to asses. 

Looking at crystallization in a system where polyte- 
trahedral frustration does not correspond to a symmet- 
ric closed-shell structure like the icosahedron would help. 
Precisely such an example is provided by the freezing of 
Ad spheres that we study in this Rapid Communication. 
It is, of course, somewhat unsatisfactory to perform a 
numerical study of a system that cannot be probed ex- 
perimentally. However, there are other examples (e.g. 
renormalization-group theory) where higher-dimensional 
model systems serve as a very useful reference state for 
the theoretical description of our 3d world. The objective 
of the numerical study that we report here is therefore 
not to present quantitative estimates of crystal nucle- 
ation barriers in Ad (even though we obtain these num- 
bers too), but to shed more light on the nature and role 
of geometrical frustration and the ease of crystallization. 

The -D4 crystal phase is formed b y st acking, without 



voids, 24-cell Platonic polytopes [16|, |17|. In general, Dd 
lattices are obtained by inserting an additional sphere 
in each void of a d-dimensional cubic lattice. In 3d the 
spacing between the spheres on the original cubic lattice 
increases to form a body-centered-cubic crystal; in Ad the 
additional sphere fits perfectly in the hole and leads to 
a unique, high symmetry crystal with maximal volume 
fraction rj = 7r 2 /16 « 0.617. There exist other dense Ad 
lattices, such as A4 and A\, but D4 packs over 10% more 
densely and offers more nearest- neighbor contacts. D 4 's 
unit cell, the 24-cell, is made of 24 octahedral cells and is 
a Platonic polytope that has no analog in other dimen- 
sions [l3| ■ Placing 24 kissing spheres around a central one 
in the 24-cell arrangement is the densest closed-shell clus- 
ter of Ad spheres [3] and is postulated to be unique [3] . 
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FIG. 1: (Color online) Equations of state of 4d hard spheres 
at constant V (closed) and P (open) extend earlier molecular 
dynamics results (2^|. Pade approximant for the fluid [2ll ] 
and Speedy fits for the crystals [22| are given for reference 
(lines). Bottom inset: at coexistence chemical potentials are 
equal, thus P CO cx = 9.15 and ^i C ocx = 13.68. Top inset: the 
common tangent to the free energy curves pinpoints the phase 
transition boundaries: ^freeze = 0.288 and ry mc it = 0.337. 



Even accounting for solvation effects, clusters with the 
24-cell geometry are locally preferred. Unlike in 3c?, for 
an equal number of particles 4c? polytetrahedral clusters 
do not form more interparticle contacts than the 24-cell, 
and their slightly larger radius offers less, not more sta- 
bilization [5[. The symmetry match between the 24-cell 
and the D4 lattice therefore guarantees that no frustra- 
tion arises from maximally kissing clusters. But neither 
the 24-cell nor any other unit cell can be assembled from 
(nearly) regular 4c? tetrahedra. Four-dimensional spheres 
are thus an ideal system to clarify the origin of geometri- 
cal frustration. An earlier compaction study of 4c? spheres 
indirectly hinted that spontaneous crystallization might 
be slow [20( , but this work could not disentangle the dif- 
ferent contributing factors, because such an analysis re- 
quires knowledge of the equilibrium phase diagram, of 
the dynamical properties of the fluid phase, and of the 
crystal nucleation barriers. Our computational study ad- 
dresses these questions. To this end, we first locate the 
4c? freezing transition, quantify the fluid order, and then 
compute the free energy barrier to nucleation at different 
supersaturations. 

Interestingly, although the equations of state of both 
the fluid and the crystal phases of 4c? hard spheres were 
computed in the early 80's [23j |. we are not aware of any 
numerical determinations of the solid-fluid coexistence 
point. Using a quasi-Maxwell construction 24j at the 
crystal stability limit [2(| , we can use these results to ap- 



proximate the coexistence range ?7coex = 0.29 — 0.35, but 
this is insufficiently accurate. To the best of our knowl- 
edge, density functional theory has only been applied to 
the fluid- A4 coexistence [25| . In order to precisely locate 
the freezing point, we thus performed standard NVT- 



FIG. 2: (Color online) Distribution of the local order correla- 
tor with I = 4 (top) and / = 6 (bottom) at P = f 9. 



Monte Carlo (MC) simulations to compute the equation 
of state of hard spheres, outside the range studied in 
Ref. [23}. As a test, we also performed constant NPT 
simulations and verified that the two techniques yielded 
consistent results. In what follows, we use the particle 
diameter a as our unit of length and the thermal energy 
ksT as our unit of energy. The equation of state of 4ci 
spheres is related to the value of the pair-distribution 
function g(r) at contact Pvo/tj = 1 + 8rig(l + ), where vq 
is the volume of a 4c? sphere and g(l + ) is the value of the 
radial distribution function at contact [2l|. The results 
for the fluid and two crystal phases are presented in Fig.[T] 
for systems containing 2048 (D4) and 4096 (fluid and A4) 
particles. The equation of state could not be calculated 
for A|, because it is mechanically unstable, which makes 
it unlikely to contribute to the crystallization process. 
We won't consider it further. To locate the fluid-solid co- 
existence regime, we need to determine the absolute free 
energy of the solid at least at one point [26[ . The absolute 
Helmholtz free energy per particle / of the D4 and A4 
crystals at ij = 0.37 is obtained by the Einstein-crystal 
method [27| • The free energy at other densities can then 
be obtained by thermodynamic integration. We find D4 
to be the thermodynamically stable crystal phase. The 
fluid- D4 coexistence pressure P COO x (Fig. Q] lower inset), 
allows to read off the melting and freezing densities by 
common tangent construction (Fig. Q] higher inset). The 
resulting two-phase region 77 coox = 0.288 — 0.337 is com- 
patible with the rough estimate above. The thermody- 
namic driving force for crystallization in the supersatu- 
rated fluid at constant pressure is the difference in chem- 
ical potential A/i = /irj 4 — /ifluid between the two phases 
displayed in the lower inset of Fig. Q] . 

To characterize the structure of the fluid and identify 
the formation of crystallites, we need a local criterion 
that distinguishes crystal from fluid. Studies in 2c? and 
3c? suggest that order parameters derived from invariant 
combinations of spherical harmonics Yi of degree / might 
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FIG. 3: (Color online) Free energy barriers for 4096 parti- 
cles at various supersaturations, along with the corresponding 
CNT and simulation parameters. The critical cluster size n* 
is obtained using a stricter order parameter (see text). 



suffice [28l.l29j. In high dimensions, it is more convenient 
to rewrite the second-order invariant in terms of Gegen- 
bauer polynomials Gf, where n = d/2 — 1, using the sum 
rule [19(. The (I + l) 2 4c? spherical harmonics give 



■ r 2 ) = yrfW E Yr(ii)Yr(h), (1) 
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where f ,■ are unit vectors. The local order correlator is 
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where the indices a and (3 run over the number of neigh- 
bors contained within a distance equal to the first min- 
imum of g(r). The local order correlation distinguishes 
between different geometrical environments: q§ set apart 
fluidlike particles from those within a D4 or an A4 lattice, 
while <74 discriminates between the two crystals (Fig. [2]). 

As freezing in 4c? is a first-order phase transition, we ex- 
pect crystallization to proceed via nucleation and growth. 
A Landau free energy analysis predicts that crystals with 
reciprocal lattice vectors forming equilateral triangles 
should initiate the nucleation [30(. Thou gh this argu- 
ment has met only limited success in 3d [29(, in Ad it 
supports the preferential nucleation of -D4, in line with 
the thermodynamic drive. To estimate the ease of crys- 
tallization, we compute the free energy barrier for crystal 
nucleation AG* . Classical nucleation theory (CNT) [3lj 
derives from the thermodynamic drive A/z and the inter- 
facial free energy 7 of a spherical crystallite a free energy 
functional that depends on the size n of the crystallite 

\{d-\)/d 



where p x is the crystal density at a given pressure and 
the shape-dependent prefactor is S4 = (1287T 2 ) 1 / 4 for 4c? 
spheres. The resulting maximal barrier height is then 
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AG*(n*) = -^—r-^ 
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AG(n) = S d (n/pjl 
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at the critical cluster size n* . The rate of nucleation per 
unit volume I is given by / = Kexp(— AG*), where k 
is a kinetic prefactor that is proportional to the diffu- 
sion coefficient in the fluid phase [29| . Though schematic 
this level of theory is sufficient to clarify the contribu- 
tion of geometrical frustration through an analysis of 7. 
Within the CNT framework the geometrical mismatch in 
3c? between icosahedral and crystal order should lead to 
a relatively large 7, while in 4c? one might expect 7 to be 
small if the locally preferred cluster scenario is valid, but 
not for polytetrahedral frustration. 

Results for 3c? crystallization are available [2^] , so only 
a few 4c? barriers are needed to complete the picture. 
Crystallization being a rare event in this regime, we 
perform constant-pressure MC runs with umbrella sam- 
pling to bias the growth of a crystal cluster from the 
fluid . A standard algorithm is employed to iden- 
tify the crystallites [H, [5^. We link pairs of nearest 
neighbors with gg > 0.4. If a particle has more than 
five links it is deemed crystalline. The number of parti- 
cles in the largest crystallite is then the order parameter. 
The resulting free energy profiles are presented in Fig. [3] 
Though qg does not discriminate between A4 and D4 
crystals, further checks with q^ show that only the latter 
nucleates. In 4c? a low q§ cutoff value is required, because 
of the minimal overlap between the crystalline and fluid 
regions (Fig. [2j) , and consequently, non-compact clusters 
are initially observed. Though the clusters irreversibly 
compactify, the process can be very slow. To reduce the 
computational burden, the system is first equilibrated by 
growing the total number of links in the largest crystal- 
lite. The low qe cutoff also artificially inflates the mea- 
sured critical cluster size. A fit to the CNT functional 
form (Eq.[3]) is thus of little use in extracting 7. However, 
because the barrier height is unaffected by this biasing 
choice and Afj, is known, we can obtain 7 from Eq. Indi- 
rectly. To validate the implied size of the CNT critical 
cluster n*-. NT we compare it to the cluster size obtained by 
imposing a purely crystalline linking criterion q^ > 0.65 
to the configurations at the top of the barrier. The dif- 
ference between the two (Fig. [3]) is no more than 25%, 
which is remarkably good in this context. 

The results of Fig. [3] allow us to conclude that the very 
slow crystallization of 4c! spheres observed in the study of 
Ref. [2fJ is due to the presence of a considerably higher 
4c? nucleation barrier than at the same supersaturation 
in 3c?. Slow nucleation could also be due to a low value of 
the kinetic prefactor k, which would require that the dif- 
fusion of particles in the dense fluid be anomalously slow. 
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But simulations with the code of Ref. [20| show no evi- 
dence for slow diffusion, not even at the highest pressures 
studied. The slow crystallization is thus consistent with a 
high degree of geometrical frustration in 4d fluids. Based 
on the similarity between the number of neighbors within 
the first peak of g(r) and the maximal kissing number 
Skoge et al. speculated that high dimensional fluids con- 
tain a number of deformed crystalline unit cells, rather 
than polytetrahedral structures [13] • However, the clear 
difference between the fluid and the 24-cell shown by the 
local order correlator (Fig [2]) suggests this not to be the 
case. The similarity between the kissing number and the 
number of first neighbors can instead be explained by a 
wide first peak of g(r) (not shown) that accommodates 
non-kissing neighbors in polytetrahedral clusters. Be- 
cause the "locally preferred" 24-cell has little to do with 
geometrical frustration, our results support the generic 
polytetrahedral structures as the source of frustration. 
By dimensional analogy, we infer that the "locally pre- 
ferred" icosahedron is not singular, but instead one of the 
many possible geometrically frustrating structures, and 
explains its limited presence in fluids. The dimension- 
less surface free-energy density r ya d ~ 1 / (fe^T) is at least 
two to three times larger in 4c? than in 3d, which indi- 
cates that geometrical frustration is surprisingly rather 
weak in 3d. It is this weakness that helps make hard 
sphere crystallization so prevalent. The interesting puz- 
zle is therefore not to identify the origin of 3d frustra- 
tion, but the source of its mildness. One possibility is 
that the tetrahedra that are part of the face-centered 
cubic (fee) structure (none are found in D4) relax the 
geometrical frustration and therefore reduce the interfa- 
cial tension. Another possibility is that the "planetary 
perturbations" that allow to exchange the positions of 
spheres at the surface of an icosahedron by sliding, go 
through a cuboctahedron configuration, which is the fee 
unit cell fl7| . If common, this phenomenon would im- 
ply that not all polytetrahedral structures are equally 
frustrating and that icosahedra might in fact be early 
nucleation sites. 

The large values for the height of the nucleation bar- 
rier of 4d crystals, as well as the evidence (Fig. [2]) that 
the local structures in the fluid and the D4 crystal are 
rather different, indicate that it is the Delaunay packing 
that matters. This finding underlines that one should be 
rather careful in caricaturing the nature of frustration as 
icosahedral in 3d liquids. Icosahedra are but one of the 
many possible polytetrahedral arrangements and little in- 
dicates that it plays a more prominent role in geometrical 
frustration than others. Note that the difficulty to crys- 
tallize 4d spheres makes them, as well as their higher 
dimensional equivalents, promising testing grounds for 
theories of packing and glass- forming liquids. 
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